HMM in Marine Sciences: Fitting a hidden Markov model to Accelerometer data
Marco Gallegos Herrada and Vinky Wang
1 Tutorial objectives
The aim of this tutorial is to delve into the process of fitting hidden Markov models (HMMs) to accelerometer data, integrating covariates into transition probabilities, and employing a classic model selection criterion to choose the most suitable model from a set of candidates. For all of this procedures, we will be using the R package momentuHMM, introduced in first part of this workshop series.
The primary learning objectives in this first tutorial are to:
- Fit a basic HMM to accelerometer data using
momentuHMM - Incorporating and interpreting covariates on behaviour transition probabilities
- Visualizing the Depth time series with decoded states
- Use the Akaike Information Criteria for model selection
- Practice what you’ve learned by working on exercises
2 Accelerometer data
Accelerometer devices measure up to three axes, which can be described relative to the body of the animal: longitudinal (surge), lateral (sway) and dorsoventral (heave). These devices are becoming more prevalent in the fields of animal biologging data as they provide a means of measuring activity in a meaningful and quantitative way. From tri-axial acceleration data, we can also derive several measures that summarize effort or exertion and relate acceleration to activity levels such as overall dynamic body acceleration (ODBA) and vectorial dynamic body acceleration (VeDBA).
ODBA and VeDBA can be used to reduce the dimensionality of three-dimension acceleration data while retaining important information. Further, because acceleration data is often at high temporal resolutions over time, it also naturally exhibits a large degree of autocorrelation, making it impossible to assume independence between sequential observations. HMMs can account for the autocorrelation present in the data while assuming that the data were generated according to a finite set of (unobserved) behaviors making them a good candidate model for this type of data structure.
2.1 Blacktip Reef Shark data
In this tutorial, we will be analyzing four days’ worth of acceleration data from a free-ranging blacktip reef shark in the Palmyra Atoll located in the central Pacific Ocean. The acceleration data was collected from a 117-cm female shark (Carcharhinus melanopterus) species using a multisensor package. This package was attached to the shark’s dorsal fin and recorded three-dimensional acceleration data at a rate of 20 Hz. It also recorded depth and water temperature at a rate of 1 Hz. After four days, the package, which was equipped with a VHF transmitter, detached from the shark and could be retrieved from the water surface (Papastamatiou et al. 2015). To assess the shark’s active behavior, the authors calculated the average ODBA over 1-second intervals, resulting in a dataset comprising 320,214 observations. Consequently, the variables in the dataset consist of time of day, water temperature (in Celsius), depth (in meters), and ODBA.
# R packages that will be used for this tutorial
library(readr)
library(momentuHMM)
library(ggplot2)
library(dplyr)
library(lubridate)
# Load data
BlacktipB <- read_delim("BlacktipB.txt",
delim = "\t", escape_double = FALSE,
trim_ws = TRUE)## # A tibble: 6 × 4
## Time Temp Depth ODBA
## <chr> <dbl> <dbl> <dbl>
## 1 7/10/2013 14:46 30.1 0.3 0.0672
## 2 7/10/2013 14:46 30.1 0.3 0.0298
## 3 7/10/2013 14:46 30.1 0.3 0.0422
## 4 7/10/2013 14:46 30.1 0.3 0.0748
## 5 7/10/2013 14:46 30.1 0.3 0.0442
## 6 7/10/2013 14:47 30.1 0.3 0.0995Exploring a quick exploration of the data, we observe there are no missing values, so we don’t have to worry about data imputation or data gaps.
## [1] 02.2 Data processing
Looking at the ODBA values throughout the observed period, we find ODBA is unusually high at some times – for this shark we assumed that values between 0 and 2 were consistent with what we expected. The dashed blue line in the plot below corresponds to a ODBA value of 2, so those values that are above this line correspond to these high ODBA values throught the observed time period:
BlacktipB_aux = BlacktipB %>%
mutate(Time = as.POSIXct(Time,format = "%m/%d/%Y %H:%M"))
BlacktipB_aux %>%
ggplot(aes(Time,ODBA)) +
geom_line() + geom_hline(yintercept = 2,linetype="dashed",color="blue")
Because accommodating extreme values can pose a problem for identification of an adequate state-dependent distribution in our HMM, we removed them from the data set. However, note that in general, deciding whether to remove extreme values or not will more likely depend on whether we find appropriate distributional forms that can accommodate them. Generally, we need to make sure that extreme values are in fact some artefact of the data collection process, not representative of a behavior of interest, or inconsistent with what we are trying to capture as well. Removing data is not good general practice but instead we can assess on a case-by-case basis.
While high-resolution data is often preferred for its ability to provide extensive information through detailed observations, it’s important to consider the potential drawbacks, as highlighted in the animal movement tutorial. One significant disadvantage of fine-scale resolution is the high autocorrelation among observations made at smaller time intervals. This increased autocorrelation can present challenges when fitting a HMM because the model may not adequately account for it. In this tutorial, we will specifically work with 1-minute average depth, temperature, and ODBA values. This significantly reduces both the time required for the model fitting process and the autocorrelation between observations. Furthermore, since we are using 1-minute averages, irregular emissions are not a concern. Finally, we will modify the variable associated with the time of day in order to be properly formatted.
# Transform into proper time format and take 1-min avg
BlacktipB_1min = BlacktipB %>% filter(ODBA <= 2) %>%
mutate(Time = as.POSIXct(Time,format = "%m/%d/%Y %H:%M")) %>%
group_by(Time = floor_date(Time, unit = "1 min")) %>%
summarise(ODBA_avg = mean(ODBA),
temp_avg = mean(Temp),
depth_avg = mean(Depth)) %>%
ungroup() %>%
mutate(Time = as.POSIXct(Time,format = "%Y-%m-%d %H:%M:%S")) %>%
as.data.frame()
head(BlacktipB_1min)## Time ODBA_avg temp_avg depth_avg
## 1 2013-07-10 14:46:00 0.05165654 30.10000 0.3000000
## 2 2013-07-10 14:47:00 0.05911751 30.10000 0.2983333
## 3 2013-07-10 14:48:00 0.04261436 30.20339 0.3271186
## 4 2013-07-10 14:49:00 0.04132418 30.81000 0.3550000
## 5 2013-07-10 14:50:00 0.04193025 31.86000 0.4050000
## 6 2013-07-10 14:51:00 0.03906949 31.31833 0.4283333After transforming the data, we now observe the 1-minute average ODBA time series across the four days in the plot below. The horizontal red lines provide a visual reference for what might represent the mean of two potential behaviors. These behaviors could be attributed to “low-activity” (for the lower values) and “high-activity” (associated with relatively higher values).Therefore, we can start proposing the model fitting with two behavioural states. It’s important to note that there are no occurrences at zero, as the purple horizontal line corresponds to the y-intercept where \(y=0\).
BlacktipB_1min %>%
ggplot(aes(Time,ODBA_avg)) +
geom_line() + geom_hline(yintercept = .05, color = "red") +
geom_hline(yintercept = .09, color = "red") +
geom_hline(yintercept = 0, color = "purple")
Now, according to the histogram of the observations associated to the ODBA values, the distribution suggest that a gamma distribution could be appropiate to model the data stream. This is reforced by the fact that we don’t have zero occurrences in this data stream, as can be visualized in the ODBA time series plot.
After selecting a data resolution and selecting a number of behaviour states, we are ready to start proposing candidate models for this data!
3 Model fitting
The best way to start when fitting a hidden Markov model is to moving from simple models to more complex ones. It’s important to remember that when fitting a hidden Markov model, we must specify two key components:
- Number of behavioural (hidden) states
- State-dependent distributions
From the latter, it’s important to to emphasize a significant assumption being made: the contemporaneous conditional independence assumption. The assumption states that, if \(\{ X_t \}\) is the observed process and \(\{ S_t \}\) is the behavioural process, for every time \(t\), the data streams are independent conditional on the hidden states. This can be mathematically expressed as follows:
\[ f(\boldsymbol{x_{t}} \mid S_{t} = s_{t}) = \prod_{p=1}^P f(x_{tp} \mid S_{t} = s_{t}),\]
where \(\boldsymbol{x_t} = (x_{t1},x_{t2},\ldots, x_{tP})\) represents the set of observations associated with the data streams at time \(t\). For the current modeling procedure, we’re dealing with a special case, since we only have one data stream. In this case, the contemporaneuous conditional independence assumption is not necessary. However, in other scenarios where the observation process includes two or more data streams, which is often the case (for instance, step length and turning angle), this assumption becomes crucial, as it enables the modeling of each data stream.
Finally, for this model fitting stage, no additonal data will be incorporated in the state-switching dynamics – which translates in not including covariances in the transition probabilities allocated in the transition probability matrix. As part of the data preparation, we need to process the data using the prepData function, which will assign the class momentuHMMData to our dataframe. This will allow us to use the functions provided by momentuHMM. Since we’re not using coordinate variables to transform them into step lengths, we need to seet the parameter coordNames to NULL.
The function fitHMM allow us to fit the model we have in mind. To perform model fitting as intended, we need to provide several arguments. The argument nbStates is used to specify the number of states we are assuming. In this case, we are working with two states. The argument stateNames is optional; it allows us to assign names to the behavioral states. In this case, we are going to name the hidden states as “low-activity” (state 1) and “high-activity” (state 2). The dist argument is used to indicate how we are modeling each of the observed data streams. In our case, this is only the ODBA values. We will model this as gamma distribution, following the earlier EDA. Moreover, when dealing with a gamma distribution for modeling observations, the standard parametrization typically involves the shape (usually denoted as \(\alpha\)) and rate (usually denoted as \(\beta\)) parameters. However, fitHMM employs the mean and standard deviation parameters instead of the shape and rate. For the choice of initial parameter values, we can use the insights gained from the previous EDA. Based on the plots presented above, setting mean values of 5 and 9 for low-activity and high-activity behaviors, respectively, for our state-dependent distributions seems to be a reasonable choice. We store these assigned initial values in the vectors mu0 for the mean and sigma0 for the standard deviation. Finally, the par0 argument requires a list of the initial values for the state-dependent distribution parameters as input.
stateNames = c("low-activity","high-activity") # define names for each state
mu0 = c(.05,.09) # initial values for the mean of each behavioural state
sigma0 = c(.02,.02) # initial values for the standard deviation of each behavioural state
fit1 = fitHMM(BlacktipBData,
nbStates=2,
stateNames = stateNames,
dist=list(ODBA_avg="gamma"),
Par0 = list(ODBA_avg=c(mu0,sigma0)))
fit1## Value of the maximum log-likelihood: 14462.08
##
##
## ODBA_avg parameters:
## --------------------
## low-activity high-activity
## mean 0.05387998 0.08484691
## sd 0.01008312 0.02486110
##
## Regression coeffs for the transition probabilities:
## ---------------------------------------------------
## 1 -> 2 2 -> 1
## (Intercept) -2.955512 -2.603854
##
## Transition probability matrix:
## ------------------------------
## low-activity high-activity
## low-activity 0.95052334 0.04947666
## high-activity 0.06889083 0.93110917
##
## Initial distribution:
## ---------------------
## low-activity high-activity
## 0.9998490135 0.0001509865Note that the value of the maximum log-likelihood is positive, which may seem unusual, as maximum likelihood values are typically negative. This positive value suggests that the likelihood surface contains very prominent peaks. Maybe data transformation can provide some insight about what’s going on with this inusual log-likelihood value. We can make use of the following property of the Gamma distribution: if \(X \sim Gamma(\alpha,\beta)\), then \(cX \sim Gamma(\alpha,\beta/c)\). Consequently, by scaling our observations up to a constant factor, the estimated parameters should correspond to those of the original data adjusted for the scaling factor. Therefore, if we attempt to fit the same model with a data transformation for the ODBA values, we should obtain identical parameter estimations, except for the scaled parameters. In this case, we chose \(c=100\) as the scale factor.
c=100
BlacktipB_1min$cODBA_avg = c*BlacktipB_1min$ODBA_avg
BlacktipBData = prepData(BlacktipB_1min,coordNames = NULL)
cmu0 = c*mu0 # initial values for the mean of each behavioural state
csigma0 = c*sigma0 # initial values for the standard deviation of each behavioural state
cfit1 = fitHMM(BlacktipBData,
nbStates=2,
stateNames = stateNames,
dist=list(cODBA_avg="gamma"),
Par0 = list(cODBA_avg=c(cmu0,csigma0)))
cfit1## Value of the maximum log-likelihood: -10244.66
##
##
## cODBA_avg parameters:
## ---------------------
## low-activity high-activity
## mean 5.388008 8.484708
## sd 1.008319 2.486117
##
## Regression coeffs for the transition probabilities:
## ---------------------------------------------------
## 1 -> 2 2 -> 1
## (Intercept) -2.955524 -2.603854
##
## Transition probability matrix:
## ------------------------------
## low-activity high-activity
## low-activity 0.95052390 0.0494761
## high-activity 0.06889082 0.9311092
##
## Initial distribution:
## ---------------------
## low-activity high-activity
## 0.9998586278 0.0001413722As expected, the estimates for the mean and the standard deviaton are the same (up to the scale factor), and all the other estimates, transtition probabilities and initial state distribution remain the same! However, note that the maximum log-likelihood is different. In this case, transforming the observed values results in ‘flattening’ the likelihood surface, but global and local maximas remain the same, therefore we have the same parameter estimates for the transition probabilities and initial state distribution. For the rest of this tutorial, we will be modelling the scaled 1-minute average ODBA values.
## $original_ODBA
## [1] -14462.08
##
## $scaled_ODBA
## [1] 10244.66# Parameter estimates for state-dependent distribution related to the ODBA data stream
list(original_ODBA = fit1$mle$ODBA_avg, scaled_ODBA = cfit1$mle$cODBA_avg)## $original_ODBA
## low-activity high-activity
## mean 0.05387998 0.08484691
## sd 0.01008312 0.02486110
##
## $scaled_ODBA
## low-activity high-activity
## mean 5.388008 8.484708
## sd 1.008319 2.486117# Parameter estimates the initial state distribution
list(original_ODBA = fit1$mle$delta, scaled_ODBA = cfit1$mle$delta)## $original_ODBA
## low-activity high-activity
## ID:Animal1 0.999849 0.0001509865
##
## $scaled_ODBA
## low-activity high-activity
## ID:Animal1 0.9998586 0.0001413722# Parameter estimates for the transtition probabilities
list(original_ODBA = fit1$mle$gamma, scaled_ODBA = cfit1$mle$gamma)## $original_ODBA
## low-activity high-activity
## low-activity 0.95052334 0.04947666
## high-activity 0.06889083 0.93110917
##
## $scaled_ODBA
## low-activity high-activity
## low-activity 0.95052390 0.0494761
## high-activity 0.06889082 0.9311092Using the plot flunction, we can visualize the state-dependent distributions using the estimated parameters. The plot below provides a depiction of low-activity and high-activity behaviors. In this visualization, the low-activity distribution exhibits higher density closer to zero, while the high-activity behavior shifts its probability mass towards higher (scaled) 1-minute average ODBA values.
## Decoding state sequence... DONE
Let’s examine the pseudo-residuals. By using the plotPR function, we can visualize these pseudo-residuals. This function displays the pseudo-residual time series, qq-plots, and autocorrelation functions (ACF). It’s important to note that pseudo-residuals are a specific type of model residuals that consider the interdependence of observations. Pseudo-residuals are computed for each data stream, and in our case, we are only considering ODBA values, resulting in a single plot of pseudo-residuals. It’s worth remembering that when the fitted model is appropriate, the pseudo-residuals should adhere to a standard normal distribution. However, based on the following plots, it appears that the model doesn’t fit well for high (scaled) ODBA values. Furthermore, the ACF plot indicates a significant autocorrelation between observations in the (scaled) ODBA time series, suggesting that the model fails to capture this autocorrelation. Notably, there is a substantial level of autocorrelation and some deviation from normality.
## Computing pseudo-residuals...
We can also compute the most probable sequence of states. As previously mentioned, the parameter estimates for transition probabilities and initial state distribution remain consistent for both the original ODBA stream and the scaled ODBA stream, resulting in identical sequences of states.
# identify most likely state using the Viterbi algorithm
BlacktipB_1min$state <- factor(viterbi(fit1))
# proportion of the behaviour states during the observed period
table(BlacktipB_1min$state)/length(BlacktipB_1min$state)##
## 1 2
## 0.5882572 0.4117428BlacktipB_1min %>% mutate(day = day(Time)) %>%
ggplot(aes(Time,ODBA_avg)) +
#facet_wrap(~day,scales = "free_x") +
geom_line(alpha=.1) +
geom_point(aes(shape=state,color=state)) + ylab("ODBA (1-min average)")
# identify most likely state using the Viterbi algorithm (scaled ODBA)
BlacktipB_1min$cstate <- factor(viterbi(cfit1))
# proportion of the behaviour states during the observed period
table(BlacktipB_1min$cstate)/length(BlacktipB_1min$cstate)##
## 1 2
## 0.5882572 0.4117428BlacktipB_1min %>% mutate(day = day(Time)) %>%
ggplot(aes(Time,cODBA_avg)) +
#facet_wrap(~day,scales = "free_x") +
geom_line(alpha=.1) +
geom_point(aes(shape=cstate,color=cstate)) + ylab("10*ODBA (1-min average)")
An essential stage in model fitting is to verify if the parameter estimates correspond to the actual Maximum Likelihood Estimation (MLE). HMMs can be susceptible to initial values, so we need a procedure to assess the reliability of the parameter estimates. The argument retryFits in fitHMM specifies the number of times initial values are perturbed before model fitting takes place. Here, we attempt to refit the model 10 times to assess whether we have identified the global maximum. This can take time to run.
set.seed(147)
fit1_s2 <- fitHMM(BlacktipBData,
nbState = 2,
stateNames = stateNames,
dist=list(cODBA_avg="gamma"),
Par0 = list(cODBA_avg=c(cmu0,csigma0)),
retryFits=10)## Attempting to improve fit using random perturbation. Press 'esc' to force exit from 'fitHMM'
##
Attempt 1 of 10 -- current log-likelihood value: -10244.66
Attempt 2 of 10 -- current log-likelihood value: -10244.66
Attempt 3 of 10 -- current log-likelihood value: -10244.66
Attempt 4 of 10 -- current log-likelihood value: -10244.66
Attempt 5 of 10 -- current log-likelihood value: -10244.66
Attempt 6 of 10 -- current log-likelihood value: -10244.66
Attempt 7 of 10 -- current log-likelihood value: -10244.66
Attempt 8 of 10 -- current log-likelihood value: -10244.66
Attempt 9 of 10 -- current log-likelihood value: -10244.66
Attempt 10 of 10 -- current log-likelihood value: -10244.66## Value of the maximum log-likelihood: -10244.66
##
##
## cODBA_avg parameters:
## ---------------------
## low-activity high-activity
## mean 5.388009 8.484710
## sd 1.008320 2.486119
##
## Regression coeffs for the transition probabilities:
## ---------------------------------------------------
## 1 -> 2 2 -> 1
## (Intercept) -2.955528 -2.603855
##
## Transition probability matrix:
## ------------------------------
## low-activity high-activity
## low-activity 0.95052409 0.04947591
## high-activity 0.06889071 0.93110929
##
## Initial distribution:
## ---------------------
## low-activity high-activity
## 9.999996e-01 3.557274e-07Seems nothing changed at all!
Now, we try a different set of initial values such that these are close to each other. Specifically, for the mean parameters, we provide values that are closer. Even when doing this, the outcome is consistent with the same model fitting, signifying that our inference process is robust.
mu2 = c*c(.04,.06) # This gives 4 and 6 as initial values for estimating the mean parameters
sigma2 = c*c(.02,.02) # This gives 2 and 2 as initial values for estimating the mean parameters
fit1_s3 <- fitHMM(BlacktipBData,
nbState = 2,
stateNames = stateNames,
dist=list(cODBA_avg="gamma"),
Par0 = list(cODBA_avg=c(mu2,sigma2)))
fit1_s3## Value of the maximum log-likelihood: -10244.66
##
##
## cODBA_avg parameters:
## ---------------------
## low-activity high-activity
## mean 5.388009 8.484710
## sd 1.008320 2.486119
##
## Regression coeffs for the transition probabilities:
## ---------------------------------------------------
## 1 -> 2 2 -> 1
## (Intercept) -2.955528 -2.60386
##
## Transition probability matrix:
## ------------------------------
## low-activity high-activity
## low-activity 0.95052409 0.04947591
## high-activity 0.06889038 0.93110962
##
## Initial distribution:
## ---------------------
## low-activity high-activity
## 9.999894e-01 1.058878e-05Finally, let’s go further by introducing a high perturbation to the initial values. Specifically, instead of 5 and 9, we will use 0.1 and 100 for the mean parameters. Additionally, for the initial values related to the standard deviation, we will set them at 100 for each state-dependent distribution. In this case, the parameter estimates remain unchanged except for one significant change: the states have swapped. In other words, the parameter initially associated with state 1 has now shifted to state 2, and vice versa. Consequently, everything has transitioned to the location of the other state, while otherwise remaining unaltered.
mu1 = c*c(.001,1)
sigma1 = c*c(1,1)
fit1_s2_long <- fitHMM(BlacktipBData,
nbState = 2,
stateNames = stateNames,
dist=list(cODBA_avg="gamma"),
Par0 = list(cODBA_avg=c(mu1,sigma1)))
fit1_s2_long## Value of the maximum log-likelihood: -10244.66
##
##
## cODBA_avg parameters:
## ---------------------
## low-activity high-activity
## mean 8.484708 5.388008
## sd 2.486117 1.008319
##
## Regression coeffs for the transition probabilities:
## ---------------------------------------------------
## 1 -> 2 2 -> 1
## (Intercept) -2.603855 -2.955526
##
## Transition probability matrix:
## ------------------------------
## low-activity high-activity
## low-activity 0.93110928 0.06889072
## high-activity 0.04947599 0.95052401
##
## Initial distribution:
## ---------------------
## low-activity high-activity
## 9.466069e-09 1.000000e+00By employing the plotPR function, we can observe that the pseudo-residuals remained consistent, despite the state parameter switching.
## Computing pseudo-residuals...
When calculating the most likely sequence of states (Viterbi sequence), it’s noticeable that the values have been reversed. This is due to the model fitting associating state 1 with “high-activity” behavior and state 2 with “low-activity” behavior.
# identify most likely state using the Viterbi algorithm
BlacktipB_1min$state_wildPar0 <- factor(viterbi(fit1_s2_long))
# proportion of the behaviour states during the observed period
table(BlacktipB_1min$state_wildPar0)/length(BlacktipB_1min$state_wildPar0)##
## 1 2
## 0.4117428 0.5882572# BlacktipB_1min %>% mutate(day = day(Time)) %>%
# ggplot(aes(Time,state_wildPar0)) + facet_wrap(~day) + geom_point()
BlacktipB_1min %>% mutate(day = day(Time)) %>%
ggplot(aes(Time,ODBA_avg)) +
#facet_wrap(~day,scales = "free_x") +
geom_line(alpha=.1) +
geom_point(aes(shape=state_wildPar0,color=state_wildPar0))
During the model fitting process, we pointed out the absence of covariates in the transition probabilities. Additionally, the ACF plot indicated a high autocorrelation between observations that the model failed to capture. Perhaps by considering the inclusion of other data streams, we could potentially reduce this autocorrelation. To explore this, we can investigate the incorporation of covariates into the state-switching dynamics to observe the potential outcomes.
4 Incorporating covariates
We can add covariates in the HMM to the transition probabilities (i.e. state process) or the state-dependent distribution (i.e. the observation process).
For example, we might be interested in including the time-of-day for each recorded observation as a covariate. We let the time of day fluctuate based on trigonometric functions, \(cos(2\pi (t/60)/24)\) and \(sin(2\pi (t/60)/24\). We construct this using the cosinor function and set period = 1440 corresponding to the minute-resolution (i.e. \(1440 = 60 \times 24\)) .
# Add tod covariate
BlacktipB_1min$tod = hour(BlacktipB_1min$Time) + minute(BlacktipB_1min$Time) / 60
tod <- BlacktipB_1min$tod
# Prep data
BlacktipBData = prepData(BlacktipB_1min,coordNames = NULL, covNames = "tod")
formula = ~ cosinor(tod, period = 1440)If our aim is to account for the effect that the time of day may have on the probability of switching between states then we incorporate the tod covariate in the transition probabilities. We demonstrate this based on the previously fitted model fit1. We use getPar0 to specify our initial parameters based on the estimated parameters of fit1.
We fit the model where the formula defines the model formula for the transition probabilities.
fit4 = fitHMM(BlacktipBData,
nbStates=2,
stateNames = stateNames,
dist=list(ODBA_avg="gamma"),
beta0 = Par0_fit4$beta,
Par0 = Par0_fit4$Par,
formula = formula)
fit4## Value of the maximum log-likelihood: 14490.46
##
##
## ODBA_avg parameters:
## --------------------
## low-activity high-activity
## mean 0.05388186 0.08506216
## sd 0.01004357 0.02485857
##
## Regression coeffs for the transition probabilities:
## ---------------------------------------------------
## 1 -> 2 2 -> 1
## (Intercept) 32.426899 -1809.14450
## cosinorCos(tod, period = 1440) -35.147711 1805.70335
## cosinorSin(tod, period = 1440) -4.270803 82.85542
##
## Transition probability matrix (based on mean covariate values):
## ---------------------------------------------------------------
## low-activity high-activity
## low-activity 0.9476722 0.05232782
## high-activity 0.1706461 0.82935388
##
## Initial distribution:
## ---------------------
## low-activity high-activity
## 1 0We can visualise how likely the shark is to be in each state at each time of day.
We observe that the shark is more likely to be in the low activity state near midday and more likely to be in the high activity state outside these hours.
Alternatively, our aim may be to account for the effect that the time of day may have on observed acceleration (i.e., ODBA) then we incorporate the tod covariate in the state-dependent distribution.
For example, we can model the mean parameter of the state-dependent Gamma distribution as a function of tod.
Building off from our fit1 model again, we use getPar0 to specify our initial parameters based on the estimated parameters of fit1.
# Fit model
## Same parameter values from fit1
Par0_fit5 <- getPar0(model = fit1,
DM = DM)
fit5 <- fitHMM(BlacktipBData,
nbState = 2,
dist=list(ODBA_avg="gamma"),
Par0 = Par0_fit5$Par,
DM = DM)
fit5## Value of the maximum log-likelihood: 14486.14
##
##
## Regression coeffs for ODBA_avg parameters:
## ------------------------------------------
## mean_1:(Intercept) mean_1:cosinorCos(tod, period = 1440)
## [1,] -3.634963 0.6941513
## mean_1:cosinorSin(tod, period = 1440) mean_2:(Intercept)
## [1,] 0.2845602 46.46127
## mean_2:cosinorCos(tod, period = 1440)
## [1,] -48.96751
## mean_2:cosinorSin(tod, period = 1440) sd_1:(Intercept) sd_2:(Intercept)
## [1,] -1.391291 -4.621587 -3.716146
##
## ODBA_avg parameters (based on mean covariate values):
## -----------------------------------------------------
## state 1 state 2
## mean 0.053566177 0.08112312
## sd 0.009837168 0.02432754
##
## Regression coeffs for the transition probabilities:
## ---------------------------------------------------
## 1 -> 2 2 -> 1
## (Intercept) -2.878229 -2.615212
##
## Transition probability matrix:
## ------------------------------
## state 1 state 2
## state 1 0.94675966 0.05324034
## state 2 0.06816577 0.93183423
##
## Initial distribution:
## ---------------------
## state 1 state 2
## 1.000000e+00 3.354248e-144We can visualise the estimated mean parameter over the time of day.
## Decoding state sequence... DONE
We see that the mean of the 1-minute averaged ODBA values vary based on the time of day for the high-activity state (state 2). In particular, the ODBA values increase from around 8am to 12am then decrease from 12-8am. The mean of the 1-minute averaged ODBA values seem to vary less with the low activity state (state 1) as it persistently increases over a 24-h period.
4.1 Model Comparison
## Computing pseudo-residuals...
## Computing pseudo-residuals...
## Computing pseudo-residuals...
There doesn’t seem to be major differences in the qq-plot and acf-plot of the pseudo-residual although there is still the residual autocorrelation and deviations in the upper tail. Perhaps including more informative covariates or increasing the number of states might help.
5 Model section: Akaike Information Criteria (AIC)
The Akaike Information Criteria (AIC) provides a measure of the trade-off between model fit and complexity. It is commonly used as a model selection criterion to select for the most suitable model among a set of candidate models. Lower AIC scores indicate a better balance between model fit and complexity.
We can compare the AIC scores using the AIC function for our original HMM (fit1), the added tod covariate in the state-dependent means HMM (fit4), and the added tod covariate in the transition probabilities HMM (fit5).
## AIC.fit1. AIC.fit4. AIC.fit5.
## 1 -28910.16 -28958.91 -28950.27We see that incorporating covariates in either the state-dependent distribution or transition probabilities lowers the AIC. fit4 has the lowest AIC which suggests that it might most suited for modelling this data among the models under consideration.
6 Exercises
- Reproduce the analysis from the blacktip reef shark data, but now consider three beahviour states instead of only two. For this, do an EDA to try to set initial values for the state-dependent distribution parameters. Look at the pseudo-residual values and the ACF plot. Does the autocorrelation reduces when considering another behavioural state.
- Is there a noticeable difference in the ACF compared to the initial 2-state HMM? \[ \]
- For the 2-state HMM, try to incorporate another variables as covariates (consider the 1-minute average temperature and depth variables). Fit this new model helping yourself from the code above.
- What are the AIC and the log-likelihood values? How are these different from the values related to the second model fitted. Based on the Akaike Information Criteria, which model seems to be the “best”?